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SUMMARY 

A simple upwind discretization of the highly coupled non-linear differential equations which define the 
hydrodynamic model for semiconductors is given in full detail. The hydrodynamic model is able to 
describe inertia effects which play an increasing role in different fields of opto- and microelectronics. A 
silicon n + — n — n + -structure is simulated, using the energy-balance model and the full hydrodynamic 
model. Results for stationary cases are then compared, and it is pointed out where the energy-balance 
model, which is implemented in most of today's commercial semiconductor device simulators, fails to 
describe accurately the electron dynamics. Additionally, a GaAs n + — n — ro + -structure is simulated 
in time-domain in order to illustrate the importance of inertia effects at high frequencies in modern 
submicron devices. Copyright © 2003 John Wiley & Sons, Ltd. 
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1. Introduction 

Our paper emerges from the the fact that today's submicron semiconductor devices are 
operated under high frequencies and strong electric fields. Information transmission using 
electromagnetic waves at very high frequencies will have a direct impact on how we design 
active and passive components in different fields of micro- and optoelectronics. In such cases, 
quasi-static semiconductor device models like the energy-balance model (EBM) are no longer 
adequate. Especially in GaAs and related materials used for high-speed device design, inertia 
effects play an important role since the impulse and energy relaxation times of the electron 
gas are close to the picosecond range. 

The most elaborate and practicable approach for the description of charge transport in 
semiconductors used for device simulation would be the Monte Carlo (MC) method pp. 
The advantage of this technique is a complete picture of carrier dynamics with reference 
to microscopic material parameters, e.g. effective masses and scattering parameters. But the 
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method must be still considered as very time consuming and hence not economical to be used 
by device designers. 

Besides the simplest concept which is the traditional drift-diffusion model (DDM), there is 
a much more rigorous approach to the problem, namely the so-called hydrodynamic model 
(HDM). This model makes use of electron temperature (or energy density) as additional 
quantities for the description of charge transport. Starting from the Boltzmann equation, 
Blotekjaer [5j and many others presented a derivation of such equations for the first time. The 
physical parameters used in the HDM can be obtained from theoretical considerations or MC 
simulations. A simplified version of the HDM is the so-called energy-balance model (EBM). 

In the first part of this work, we give a short definition of the charge transport models. We 
will illustrate how the different models and their parameters can be related to each other. In 
a second part, we give a simple discretization scheme for the full hydrodynamic model. In the 
last part, we compare the different models for the case of a submicron silicon ballistic diode 
(an n + — n — n + structure) and a Gallium Arsenide ballistic diode. 



2. The HDM for silicon 

Since we will illustrate the HDM for the case of an n-doped ballistic diode where the 
contribution of electron holes to the current transport is negligible, we will only discuss the 
charge transport models for electrons. Generalization of the models to the case where both 
charge carriers are present is straightforward. 

The four hydrodynamic equations for parabolic energy bands are 

| + VJ=0 , (1 > 

1 £ + (Vj)v + (jV)v = nE--V[ , 2 

at m m \ e J t p 

^ j-^ jjji 'y 'j. 

-envE - V(nkTv) - V(-kVT) ? , (3) 

and 

V(eV$) = e(n - N D ) , (4) 

where n is the electron density, v the drift velocity of the electron gas, e > the elemental 
charge, E = — V$ the electric field, $ the quasi-static electric potential, T the electron gas 
temperature, to the electron energy density, k the thermal conductivity, e is the dielectric 
constant, and Nd is the density of donors. 

The particle current density j = nv is related to the current density J by the simple formula 
J= -ej. 

Eq. is simply the continuity equation which expresses particle number conservation. Eq. 
J21 is the so-called impulse balance equation, eq. Q the energy balance equation and eq. (@J 
the well-known Poisson equation. 
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We will solve the hydrodynamic equations for n,j,u> and $. To close the set of equations, 
we relate the electron energy density to the thermal and kinetic energy of the electrons by 
assuming parabolic energy bands 

lo = -nkT H — nmv 2 . (5) 
2 2 w 

In fact, we already assumed implicitly parabolic bands for the impulse balance equation, which 
is usually given in the form 

^+{Wp)v + (pV)v = -enE-W(nkT)- , (6) 

Ot T p 

i.e. we replaced the electron impulse density by the particle current density by assuming 
p = mnv, where to is a constant effective electron mass. 

The impulse relaxation time t p describes the impulse loss of the electron gas due to the 
interaction with the crystal, the energy relaxation time t w the energy transfer between the 
electron gas with temperature T and the crystal lattice with temperature Tr,. t p and t u are 
usually modelled as a function of the total doping density Np + Na where Na is the density 
of acceptors, the lattice temperature Tl, the electron temperature T or alternatively the mean 
energy per electron u)/n. 

A simplification of the full HDM is the energy-balance model. In the EBM, the convective 
terms of the impulse balance equation are skipped. The energy balance equation is simplified by 
the assumption that the time derivative of the mean electron energy du/dt is small compared 
to the other terms and that the kinetic part in lo can also be neglected, i.e. 

lo = -nkT . (7) 
2 v ; 

This non-degenerate approximation which avoids a description by Fermi integrals is justified 
for the low electron densities in the relevant region of the simulation examples, where velocity 
overshoot can be observed. 

The energy balance equation then becomes 

V(^nfcTw) - V(kVT) = 

^ §nfc(T-T L ) 
-envE - ^ — — , (8) 



e ^t/nkl \ 

— T p nE r p V . 9 

to V e / 



and the impulse balance equation becomes the current equation 

e -± e ^(nkT^ 
— T„nE 1 

TO TO 

Continuity equation and Poisson equation are of course still valid in the EBM. Neglecting 
the time derivative of the current density is equivalent to the assumption that the electron 
momentum is able to adjust itself to a change in the electric field within a very short time. 
While this assumption is justified for relatively long-gated field effect transistors, it needs to 
be investigated for short-gate cases. 

A further simplification of the EBM leads to the drift-diffusion model. The energy balance 
equation is completely removed from the set of equations, therefore it is no longer possible to 
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include the electron temperature T in the current equation. T is simply replaced by the lattice 
temperature Tl- Therefore the DDM consists of the continuity equation, the Poisson equation 
and the current equation 



J 



CT, 



^nE 



CT, 



. Vn 

m V e / 



(10) 



and it is assumed that the electron mobility is a function of the electric field. But at least, the 
electron temperature is taken into account in an implicit way: If one considers the stationary 
and homogeneous case in the HDM, where spatial and temporal derivatives can be neglected, 
one has for the current equation 



J 



-nE 



or 



m 

6 -^E 
m 



and the energy balance equation becomes simply 



^ |fc(T - T L ) 
—evE — — 



1 9 

pi) 



(11) 

(12) 
(13) 



-k(T — T L ) 



(14) 



Combining eq. (|12|) and i|13|) leads to the relation 

^(T p T^-T p 2 /2)£ 2 = fi 

m F 2 

In our simulations for silicon, we will use the Baccarani-Wordemann model, which defines the 
relaxation times by 

t p = —fJ-ajr ( 15 ) 



m 
Ye 1 



3 k 



T 
1 

2 + 2^?T 



2 eVg 
3fc 



TT 

t + t 

T 2 



Ti 



(16) 



v s is the saturation velocity, i.e. the drift velocity of the electron gas at high electric fields. 
jiQ is the low field mobility, which depends mainly on the lattice temperature and the total 
doping density. 

For the sake of completeness, we mention that inserting the expressions for the relaxation 
times into eqs. l(T2*)l and (|T3)l leads to the 2?(T)-relation 



E 2 = \ 
Mo 



\T L ) 



and the electron mobility [i — (e/m)T p is given by 

er p (E) _ 



KE) 



Mo 



rn 



(17) 



(18) 



1 
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This is the well-known Caughey-Thomas mobility model It has the important property 
that v(E) ~ u Q E for £«^andi;~ v s for E » 

v ' ^ Mo Mo 

The EBM has the big advantage that it includes the electron temperature T, such that the 
electron temperature gradient can be included in the current equation, and the mobility can 
be modelled more accurately as a function of T. 

An expression is needed for the thermal conductivity of the electron gas, which stems from 
theoretical considerations 

k= (5/2 + r> fc2/i(T) T . (19) 
e 

Several different choices for r can be found in the literature, and many authors |S] even 
neglect heat conduction in their models. But Baccarani and Wordeman point out in [Sj that 
neglecting this term can lead to nonphysical results and mathematical instability. Although 
their work is directed to Si, their remarks should be equally valid for GaAs since the equations 
have the same form in both cases. We will present a GaAs MESFET simulation comparable 
to the one of Ghazaly et al. jS] in a forthcoming paper, but with heat conduction included. 
The best value for r appears to be —2.1 for silicon at 300 K, according to comparisons of 
hydrodynamic and MC simulations of the ballistic diode 0. 

3. Discretization scheme 

Today, many elaborate discretization methods are available for the DDM equations or EBM 
equations. The well-known Scharfetter-Gummel method ||] for the DDM makes use of the fact 
that the current density is a slowly varying quantity. The current equation is then solved exactly 
under the assumption of a constant current density over a discretization cell, which leads to 
an improved expression for the current density than it is given by simple central differences. It 
is therefore possible to implement physical arguments into the discretization method. Similar 
techniques have been worked out for the EBM J0| • But due to the complexity of the HDM 
equations, no satisfactory discretization methods which include physical input are available for 
this case. For the one dimensional case, this is not a very big disadvantage, since the accuracy 
of the calculations can be improved by choosing a finer grid, without rising very strongly the 
computation time. 

Therefore, we developed a shock-capturing upwind discretization method, which has the 
advantage of being simple and reliable. For our purposes, it was sufficient to use a homogeneous 
mesh and a constant time step. But the method can be generalized without any problems to 
the non-homogeneous case. Also a generalization to two dimensions causes no problems, and 
we are currently studying two dimensional simulations of GaAs MESFETs. 

The fact that the discretization scheme is fully explicit should not mislead to the presumption 
that it is of a trivial kind. In fact, stabilizing a fully explicit discretization scheme for such a 
highly non-linear system of differential equations like the HDM is a difficult task, and a slight 
change in the discretization strategy may cause instabilities. Therefore, naive application of 
the upwind method does not lead to the desired result. The order how the different quantities 
are updated is also of crucial importance for the maximal timesteps that are allowed. The 
timesteps can be enhanced by using an implicit scheme, but only at the cost of an increased 
amount of computations needed for the iterative numerical solution of the implicit nonlinear 
equations. 
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The device of length I is decomposed into N cells C$ of equal length Ax — l/N . 'Scalar' 
quantities like the electron density rii, (i = 1, ...N), the potential <I>i and the electron energy 
density Wj are thought to be located at the center of the cells, whereas 'vectorial' quantities 
like the particle current Jj+1/2, {i = 1, ..JV — 1) and the electric field E i+ i/ 2 — ($j — $ i+1 )/Aa; 
are located at the boundaries of the cells. If necessary, we can define e.g. E^ by 

E i = ^(E i _ 1/2 +E i+1/2 ) , (20) 

but a different definition will apply to e.g. ji, as we shall see. 

The fundamental variables that we will have to compute at each timestep are n;, u>i (or 
Ti) for i = 2, „JV — 1, and ji+1/2 for i = 1, „JV — 1, if m, n^r, $1, <&at, wi, and ujn are fixed 
by boundary conditions. All other variables used in the sequel should be considered as derived 
quantities. 

The constant timestep At used in our simulations was typically of the order of a few tenths 
of a femtosecond, and quantities at time T — tAt carry an upper integer index t. 

Having calculated n\ for a timestep t, we define the electron density at the midpoint by 



t _ J l n i \ n \-i '■ il+1/2 > 

i+l/2 — \ 3 t _ I t .At < 



(21) 



i.e the electron density is extrapolated from neighbouring points in the direction of the electron 
flow, and further 

V i+l/2 — ji+1/2 / n l+l/2 ■ (22) 

The upwind extrapolation of the electron density which is given by the weighting factors 3/2 
and —1/2 is improving the accuracy of the scheme compared to the usual upwind choice 
n i±i/2 — n i, where simply the neighbouring value in upwind direction is used. Analogously we 
define 

f 3 t _ 1 At \ At >n 

■t _ ) 2 J i-l/2 2 J i-Z/2> ■ Ji+1/2 ^ u / O o> 

Ji — \ 3 At 1a \ ■ At ^ n ' 



and 



2 ji+1/2 ~ 2^+3/2) : Ji+1/2 < 



vl=3lK ■ (24) 



The discretization of the Poisson equation can be done by central differences. The continuity 
equation is discretized as follows: 

l£l_A = J+^-^ i = 2 ,...N-l , (25) 
At Ax y ' 

and thus can be calculated from quantities at T = tAt. Eq. defines a conservative 
discretization, since the total number of electrons can only be changed at the boundaries, 
where electrons may enter or leave the device. Electrons inside the device which leave cell Ci 
at its right boundary enter cell Cj+i from the left. The values of Jn+i/2 an d J1/2 will not be 
needed in our simulations, since we will use boundary conditions for the electron density which 
fix n\ and tin- 
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As a next step we have to discretize the impulse balance equation. Most of the terms can 
be discretized by central differences: 



■4+1/2 -h+i/2 _ _ e , t 

At ~ m ^+V2^+l/2 

hn\ +1 T} +1 -nW)/Ax-±^-{conv)\ +1/2 , (26) 

171 T p,i+l/2 

Et +1/2 = m-^ i+1 )/Ax , (27) 
but the convective terms require an upwind discretization 

(conv)l +1/2 = jt +1/2 (vl +1/2 -vl_ 1/2 )/Ax + 

«l + i/ 2 0'i + i/2-.?ti/2)/Az (28) 

if , 2 or w* +1 / 2 have positive direction and otherwise 

{conv)\ +1/2 = jl +1/2 (vl +3/2 -vl +1/2 )/Ax + 

^+1/2(^+3/2-^+1/2)^. (29) 



We observed that the stability of the scheme is improved for silicon if the current density is 
first updated by 

;■*+! _ at 

At ~ --^+1/2^+1/2 - 

hnt +1 T! +1 -r4Tl)/Ax-^- , (30) 

and then j* +1 / 2 i s updated by the convective terms, but with j\ +l / 2 and w* +1 / 2 replaced by 
the values resulting from i* +1 / 2 - 

The electron temperature is related to the energy density by the relation wf = \nkTl + 
^■mn\v t2 and can therefore be regarded as a dependent variable. The energy balance equation 
is discretized by defining first 



1 i w *-3 w i-i : i'+i/ 2 >o 

0J i+i/2 — S 3, ,t i,,t .At ^ n ' ^ 

I 2 W i+l ~ 2 UJ *+2 ■ h+l/2 < U 



such that Tl +X j 2 is also defined by our upwind procedure. Then the discretization is given by 

= -en\v\E] 4-^ 



At 



£ x (jl,i+l/2 ie,i-l/ 2 ) 

~Ai^' i+1/2 ~~ 4,1-1/2) 

— ^j0'fc,i+l/2 ~ 3h,i- 1/2) i (32) 

Copyright © 2003 John Wiley & Sons, Ltd. tot. J. Numer. Model. 2003; 16:161-174 

Prepared using jnmauth.cls 



168 



A. ASTE AND R. VAHLDIECK 



where we have defined three energy currents 

ie,i+l/2 = v i+l/2 UJ i+l/2 ' (33) 

3p,i+l/2 — ^ii+l/2^i+l/2 ) (34) 

and 

3ii+i/2 = -< r T* +1 -Ti)/Ax . (35) 

4. Stationary simulation results 

We simulated an n + — n — n + ballistic diode, which models the electron flow in the channel of 
a MOSFET, and exhibits hot electron effects at scales on the order of a micrometer. Our diode 
begins with an 0.1 /im n + "source" region with doping density Nd = 10 18 cm -3 , is followed 
by an 0.1 /im n "channel" region (Nd = 2 • 10 15 cm~ 3 ), and ends with an 0.1 /im n + "drain" 
region (again Nd — 10 18 cm~ 3 ). The doping density was slightly smeared out at the junctions. 
We used the following physical parameters for silicon at Tl = 300K The effective electron 
mass m = 0.26TO e , where m e is the electron mass, e = 11.7, and v s — 1.03 • 10 5 m/s. The low 
field mobility is given by the empirical formula 

(M>(N D ) = /w + 1 + {N ^ Nref)0 . 72 , (36) 
fx min = 80cm 2 /Vs , Afi = 1430cm 2 /Vs - fi min , (37) 

N ref = 1.12 ■ 10 17 cm- 3 . (38) 

The temperature dependent mobilities and relaxation times follow from the low field values 
according to eqs. (| 1 51 1 6|) . 

For boundary conditions we have taken charge neutral contacts in thermal equilibrium with 
the ambient temperature at x = and x = I = 0.3 /im, with a bias V across the device: 

m = N D (0) , n N = N D (l) , (39) 

T 1= T N = T L , $ 1 = , $ N = V . (40) 

Iinitial values were taken from a simple DDM equilibrium state simulation. 

Stationary results were obtained by applying to the device in thermal equilibrium a bias 
which increased at a rate of typically 1 Volt per picosecond from zero Volts to the desired final 
value. After 6 picoseconds, the stationary state was de facto reached, i.e. the current density 
was then constant up to 10~ 4 %. 

In most cases we used N = 200 discretization cells, which proved to be accurate enough, 
and time steps At of the order of a femtosecond. A comparison with simulations with N > 500 
shows that all relevant quantities do not differ more than about 5% from the exact solution. 

The computation of a stationary state on a typical modern workstation requires only few 
seconds of CPU time, if FORTRAN 95 is used. 
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Figure 1. Electron density for V=l V (stationary state). The HDM exhibits more structure than the 
EBM. The symmetric curve is the nearly abrupt doping profile. 



Fig. 1 shows the electron density for the EB and HD charge transport models for a bias of 
1 V. The choice of 1 V is meaningful since at higher bias (> 2 V), the HDM would no longer 
be applicable or the device would even be destroyed. 

In this paper, solid lines refer always to the HDM anddashdotted lines to the EBM. The 
HDM exhibits more structure than the EBM. It is interesting to observe that the electron flow 
becomes supersonic at x = 0.109 /im in the HDM (whereas it remains subsonic in the DDM). 
Fig. 2 shows also the soundspced in the electron gas calculated from the electron temperature in 
the HDM, which is given by c = y/ kT/m if heat conduction is included and by c = \JhkT/3rn 
otherwise (dashed curve). In fact, a shock wave develops in the region where the Mach number 
v/c is greater than one. In the DDM, the electron velocity exceeds the saturation at most by 
30%. The maximum electron velocity in the HDM is 2.61 v s , in the EBM only 1.87 v s . 

Finally we observe in Fig. 3 that the EBM is able to describe the electron temperature in 
an acceptable way. It also predicts the cooling of the electron gas near x = 0.1 /*m, which is 
caused by the little energy barrier visible in Fig. 4, where the electric field has a positive value. 
The dotted line in Fig. 4 show the electric field in thermal equilibrium. 

But still it is interesting to note that the drift part of the mean electron energy u> jn becomes 
large in a small region around the drain-source junction, where the electron kinetic energy can 
be as large as 58% of the total energy (Fig. 5). 

From the point of view of device modeling, the J-V-characteristics resulting from the three 
models is of importance (Fig. 6). At low bias, the electron mobility is in all three models is 
governed by the low field mobility fi . 

But it is quite astonishing how the EBM predicts a J-V-curve which is in very good agreement 
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Figure 2. Electron velocity for the different charge transport models (V=l V). Between x = 0.109 fim 
and x = 0.146 fim the HDM electron velocity is supersonic. The dashed curve is the soundspeed. 
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Figure 3. Electron temperature for V=\ V. HDM and EBM are in good agreement. (The DDM electron 
temperature calculated from the electric field becomes meaningless in this case.) 
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Figure 4. Electric field inside the device. 
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Figure 5. Ratio of kinetic drift energy and total energy of electrons for V=l V. 
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Figure 6. J — ^-characteristics of the device (HDM). The dots represent data obtained from the EBM. 

with the HDM prediction also in the range of higher biases, such that the two curves in Fig. 
6 are nearly undistinguishable. As we have already mentioned, the EBM does not take inertia 
effects into account, which play no role for the stationary case. The predictive power of the 
two models is quite different, as we will see in the next section. But in the stationary case, the 
difference in the J-V-characteristics is, roughly speaking, averaged out. 



5. Inertia effects 

For GaAs, the relaxation times are quite high. Therefore, inertia effects will become important 
if the applied electric field changes at a high frequency. We simulated therefore a GaAs ballistic 
diode at 300 K. The diode begins with an 0.2 fim n + source region with doping density 
Ne> = 2 • 10 17 cm~ 3 , is followed by an 0.4 /im n channel region (Nn = 2 • 10 15 cm~ 3 ), and 
ends with an 0.2 ^m n + drain region (No = 2 • 10 17 cm~ 3 ). The relevant data like energy- 
dependent relaxation times and electron mass were obtained by two-valley MC simulations, 
where also the non-parabolicity of the two lowest conduction band valleys in GaAs was taken 
into account. For the sake of brevity we will not go into details here, which will be given in a 
forthcoming paper concerning the full hydrodynamic simulation of a GaAs MESFET structure. 

In order to show the different behavior of the HDM and EBM in time-domain, we applied to 
the 0.1 V pre-biased ballistic diode an additional 0.1 V pulse of 1 picosecond duration (see Fig. 
7) . Fig. 8 shows the particle current density in the exact middle of the device for both models 
as a function of time. Whereas the current in the EBM reacts immediately to the applied field, 
the current in the HDM shows relaxation effects. The impulse relaxation time in the channel 
of the diode is of the order of 0.3 picoseconds, the energy relaxation time lies between 0.3 and 
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Figure 7. Applied bias as a function of time. 



1 picosecond. We emphasize the fact that considering the total current 



Jtot 



dE 

-env + e e r — - 
at 



(41) 



does not help; the effect remains. Therefore we must conclude that the EBM, which is often 
termed " hydrodynamic model" in commercial semiconductor device simulators, may lead 
accidentally to reasonable (static) characteristics of a device, although the physical processes 
inside the device are modelled incorrectly. 



6. Conclusion 

A simple discretization scheme for the hydrodynamic model is given in full detail which 
gives a valuable tool to the practitioner entering the field of hydrodynamic device modeling. 
Comparisons of the different transport models show that the energy-balance model is capable 
of describing the behavior of submicron devices fairly well, but full hydrodynamic simulations 
are needed in order to give a satisfactory description of the device from the physical point 
of view. At high frequencies, inertia effects become important in GaAs and related materials. 
It is therefore clear that the EBM will be no longer adequate for simulation of high-speed 
submicron devices in the near future. 
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